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ABSTRACT 

We present simulations of scattering phenomena which are important in pulsar observations, but 
which arc analytically intractable. The simulation code, which has also been used for solar wind 
and atmospheric scattering problems, is available from the authors. These simulations reveal an 
unexpectedly important role of dispersion in combination with refraction. We demonstrate the effect 
of analyzing observations which are shorter than the refractive scale. We examine time-of-arrival 
fluctuations in detail: showing their correlation with intensity and dispersion measure; providing a 
heuristic model from which one can estimate their contribution to pulsar timing observations; and 
showing that much of the effect can be corrected making use of measured intensity and dispersion. 
Finally, we analyze observations of the millisecond pulsar J0437— 4715, made with the Parkes radio 
telescope, that show timing fluctuations which are correlated with intensity. We demonstrate that 
these timing fluctuations can be corrected, but we find that they are much larger than would be 
expected from scattering in a homogeneous turbulent plasma with isotropic density fluctuations. We 
do not have an explanation for these timing fluctuations. 
Subject headings: pulsars: general - ISM:general 



1. INTRODUCTION 

Observations of radio pulsars have shown the effects 
of the interstellar plasma since pulsars were discovered. 
The first and most obvious effect is dispersion due to 
the column density of free electrons between the pulsar 
and the observer. However the effects of scattering due 
to fluctuations in the electron density were soon recog- 
nized as they are very strong in pulsars (Rickett, 1969; 
Rankin et al., 1970; Cordes, 1986). It is now known 
that scattering by fluctuations in electron density due to 
Kolmogorov turbulence are diffractive in nature but the 
diffractive process is strongly modulated by large scale 
refraction. As a result there are two spatial scales of 
the intensity fluctuations, a small scale now called the 
"diffractive scale" Sdif and a larger one called the "re- 
fractive scale" s rc f (Prokhorov et al., 1975). The early 
observations showed only the diffractive scale. The re- 
fractive scale was discovered by Sieber (1982) and ex- 
plained by Rickett et al. (1984). Furthermore the col- 
umn density changes slowly and this is also observable 
(Rawley et al., 1988; Ramachandran et al., 2006; You 
et al., 2007). These fluctuations in "dispersion measure" 
(the pulsar observer's term for column density) have their 
origin in the same interstellar turbulence and they merge 
with the diffractive and refractive effects. Observations 
and theory of interstellar scattering have been reviewed 
by Rickett (1990) and Narayan (1992). 

Scattering is an inherently spatial effect, but it is ob- 
served as time variation in, for example, the pulse in- 
tensity. The temporal variation is simply caused by spa- 
tial variation convected past the observer by the relative 
motion of the pulsar, the Earth, and the ionized inter- 
stellar medium (IISM) . Thus we relate an observed time 
scale r to the corresponding spatial scale S using an ef- 
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fective velocity S = V c tfT (Cordes and Rickett, 1998). 
As there are two spatial scales Sdif and s rc f, there are 
two temporal scales Tdif and r rc f. In practice Tdif is usu- 
ally minutes to hours so the diffractive component is well 
sampled in a typical observation. However r ro f is usually 
days to months, so refractive variations are seldom seen 
in a single observation period. This is why they were 
not discovered for nearly two decades. Consequently ob- 
servations are usually analyzed neglecting the refractive 
variation, but the derived parameters of the diffractive 
effects are seen to vary with observing epoch (Gupta et 
al., 1994; Gupta et al., 1999). The question then arises, 
"Is this variation part of a continuous variation due to 
homogeneous turbulence in the IISM, or must we invoke 
an inhomogeneity in the electron density as was used by 
Fiedler et al. (1987) to explain extreme scattering events 
(ESE)?" 

2. SIMULATION OF SCATTERING 

Scattering theory only provides asymptotic solutions 
for the moments of the electric field in general and these 
are not sufficient for comparison with most observations, 
nor are they adequate to predict the effect of analyzing 
observations on a time scale which is shorter than the re- 
fractive time. Simulations have long been used to study 
the diffractive and refractive scales (Coles & Filice 1984; 
Coles et al. 1995a), but they have not, until this work, 
been extended to include the still larger scales of disper- 
sion measure variations, to study the effects of analyzing 
relatively short observations or to study time-of-arrival 
variations. The code used in this work, which is based 
on that used in the studies referenced above, is publicly 
available and may be obtained from the authors (contact 
bcoles@ucsd.edu). 

The scattering process can be simulated by decompos- 
ing the medium into a series of thin screens. Turbu- 
lence in these screens can be generated numerically and 
the incident wave propagated from screen to screen nu- 
merically. The propagation between screens can be for- 



mally written as a Fresnel diffraction integral, which is 
a two-dimensional convolution. For numerical work it is 
attractive to implement this convolution in the Fourier 
tranform domain. The electric field is transformed into 
an angular spectrum, the angular spectrum is propagated 
to the next screen by multiplying each component by the 
appropriate propagation constant, and finally the angu- 
lar spectrum is inverse transformed to provide the electric 
field incident on the next screen. 

These calculations can be repeated at different observ- 
ing frequencies with the same screen (choosing cither 
plasma dispersion or non-dispersive frequency scaling), 
thus building up a data-cube E(x, y, f) in the observing 
plane. With frequency information one can recover pulse 
shapes and simulate dynamic spectra. 

In practice it appears that scattering observed in many 
pulsars is dominated by a compact region and thus can 
be modeled by a single thin screen. In this case there 
can be no velocity shear in the medium and the effec- 
tive velocity is well defined. Here we consider only the 
canonical case of a plane wave incident on a single thin 
screen. The important case of a spherical wave incident 
on a thin screen can be derived from these results with 
the appropriate scaling, as given by Rickett et al. (2000) 
in their Appendix A. 

The details of the simulations used in this paper are 
thoroughly discussed by Coles et al. (1995a) and the 
interested reader is referred to this exposition. For our 
purposes here it is necessary to realize that the scattering 
region is represented by a finite grid. The Fourier trans- 
form, implemented discretely, requires that the grid be 
periodic. Thus we must make sure that the period or 
"window" (L x x L y ) is sufficiently large to adequately 
represent the largest scale of the process s re f, yet suffi- 
ciently finely sampled (dx, dy) to catch the smallest scale 
of the process Sdif- As the scattering strength increases 
Sdif and s re f separate further. Their geometric mean re- 
mains the "Fresnel scale" r{. Thus r r = y/z/k provides 
an appropriate normalizing scale. Here k = 2tt/X is the 
wavenumber and z is the distance from the scattering 
screen to the observer. To first order the window should 
be adjusted so r{ is the geometric mean of dx and the 
smaller of L x and L y . This separation of scales means 
that the number of samples required for the grid increases 
like the fourth power of the strength of scattering and 
this sets a rather firm bound on the strongest scattering 
that it is feasible to simulate. The nature of this limit, 
expressed in terms of the number of samples required to 
calculate the field with a given accuracy, is discussed by 
Coles et al. (1995a). Fortunately in very strong scatter- 
ing one can often use asymptotic expressions. 

In this paper we will examine several specific ques- 
tions of interest. First, we will consider how the param- 
eters that one would estimate for diffractive scintillation 
from relatively short observations, will vary over much 
longer time scales. Second, we will consider the phe- 
nomenon of "tilted structure" in dynamic spectra and 
the related asymmetrical parabolic arcs, which have been 
attributed to dispersive refraction (Ewing et al., 1970; 
Shishov, 1974; Hewish et al., 1985; Cordes et al., 2006). 
We will confirm that this is caused by a combination of 
dispersion and refraction - it is not apparent in a simula- 
tion with non-dispersive refraction. Third, we will show 



how the time-of-arrival (TOA) varies due to scattering. 
We will show the correlations between TOA variation, 
intensity and dispersion measure. We will demonstrate, 
using both simulations and observations of the pulsar 
J0437— 4715 made at Parkes, that much of the observed 
TOA variation can be corrected if measurements of in- 
tensity and dispersion measure are available. We provide 
heuristic formulae, derived from the simulations, from 
which the TOA variations for a given pulsar can be esti- 
mated. We will not, in this paper, discuss the effects of 
inhomogeneous or anisotropic turbulence although our 
simulation code is well-suited to such studies and they 
may be important for many pulsars. 

3. SCATTERING THEORY 

Scattering in the IISM is well-modeled as small-angle 
forward-scattering by irregularities in refractive index 
which are very large compared with the wavelength. In 
this case one can write parabolic partial differential equa- 
tions for the moments of the electric field under the as- 
sumption that the scattering medium is delta correlated 
in the propagation direction (z). This is sometimes called 
the Markov assumption. This theory has been known 
for some time (Tatarski 1967; Gochelashvily & Shishov; 
1971; Prokhorov et al. 1975), but solutions to the equa- 
tions have been limited. The second moments of the elec- 
tric field are relatively well-behaved, so we have a good 
expression for the brightness distribution and a reason- 
able expression for the impulse response. However the 
fourth moments, which describe the intensity statistics, 
are more difficult. They can be solved in weak scintilla- 
tion using the Born approximation, and an asymptotic 
solution is available in strong scintillation. In moderately 
strong scintillation one can develop asymptotic series ap- 
proximations but this is difficult and the results are hard 
to scale. One can use numerical simulations, or numerical 
solutions to the moment equations in this region. Good- 
man and Narayan (2006) have used the latter approach 
and gave a set of empirical approximation formulae valid 
in this transition region. Observations often reveal phe- 
nomena which are not easily expressed as moments of 
the electric field. It is difficult to describe their statis- 
tics analytically, even in weak scintillation. It is in these 
cases that simulations are presently indispensable. 

Here we will summarize the widely-used equations for 
quantities, such as the spatial scales of the field, which 
arc necessary to interpret observations. Most of the 
widely used results are for asymptotically strong scin- 
tillation and one of the objectives of this work is to see 
where these equations begin to break down in moder- 
ately strong scintillation. We will limit our discussion to 
isotropic power-law spectra of refractive index for which 
the exponent is close to the Kolmogorov value (-11/3). 
We will consider only the case of a plane wave incident on 
a thin screen of homogeneous spatial fluctuations which 
do not change with time. This is a canonical case which 
can be mapped into the case of a spherical wave incident 
on the screen, and more complex problems can be assem- 
bled using multiple screens either with a plane wave or a 
spherical wave incident. 

The screen is a thin region in which the refractive index 
differs from its mean value by Sn. This screen changes 
the phase of an incident plane wave by <j)(r) = k 5n(r) dz. 
Here k = 2tt/X is the spatial wave frequency and Sz is 



the screen thickness. In a plasma Sn oc A 2 so <p(r) oc A. 
In a non-dispersive medium Sn is independent of A so 
<fr(r) ex 1/A. The phase variations are best described, for 
our purposes, by a structure function 



m 2 = 2 Q r(l + a/2) cos{air/2)D(r f ). 



(3) 



D(s) = <(0(r)-0(r + s)) 2 > 



(1) 



The structure function exists if the random process 
4>{r) has stationary differences. This is a weaker condi- 
tion than the wide-sense stationarity which is necessary 
for the existence of a covariance function. It is particu- 
larly useful for processes with power law spectra, where 
the spectral exponent (3 is flatter than -4. In such cases 
the structure function is also power law D(s) = (s/so) a , 
with a — — (/3 + 2). Thus the structure function ap- 
plies to Kolmogorov turbulence in the inertial sub-range 
(/3 = —11/3) and it is also attractive for propagation 
calculations because it arises naturally in the partial dif- 
ferential equations for the moments of the electric field. 
The structure function exists even if the outer scale has 
not been defined and the variance diverges. 

We normalize the electric field so the mean intensity is 
unity. The transverse spatial correlation of the electric 
field at the output of the screen for a monochromatic 
plane wave incident on the screen can be written in any 
strength of scintillation as 

T E (s) = <E(r)E*(r + s)> = exp(-0.5 D(s)). (2) 

This result is independent of the distance from the 
screen. In the power law case, with < a < 2, we 
have r#(s) = exp(— 0.5 (s/s ) Q ), so the e -1 / 2 scale 
of the field is so the spatial separation at which the 
rms phase difference is one radian. In the language of 
interfcrometry Te(s) is the visibility. The brightness 
distribution B(k9) is just the spatial Fourier transform 
F{.Te(s)}. Of course B(9) is also independent of dis- 
tance from the screen (z). Both Te(s) and B(9) are 
quasi-Gaussian so the e -1 / 2 width of B(9) is 6q w 1/kso- 
In fact both functions have long tails which are often 
important, but their e -1 ' 2 widths are quite close to the 
Gaussian approximation. The impulse response I(t) can 
be derived from B(9) as I(t)dt = 2ir9B(9)d9, where 
t = 9 2 z/2c (for a plane wave incident). This yields 
I(t) = (2-kcJ z)B(9 = y / 2tc/z). In the case of a quadratic 
structure function, thus a Gaussian B(9), one obtains 
I(t) = Jo exp(— i/*o) where t = 9 2 z/c. The e _1 width 
of I(t) is close to to for Kolmogorov spectra. This is also 
valid in any strength of scintillation, requiring only that 
the screen be thin. 

The intensity covariance Tx(s) = <5T(r)5T{r + s)>, 
where 51 = X— < I > , can be calculated in weak scin- 
tillation where a closed form expression for its Fourier 
transform can be derived using the Born approximation. 
The Born approximation is valid where the intensity vari- 
ance is less than unity. We use the Born variance m 2 as 
a measure of the strength of scintillation. When it is less 
than unity the scintillation is weak and it is a good ap- 
proximation to the actual variance. When m 2 is large the 
scintillation is not weak and it is not a good approxima- 
tion to the actual variance, but it remains a very useful 
measure of the amount by which the scintillation exceeds 
the weak condition, i.e. the strength of scintillation. The 
intensity spectrum in weak scattering can be written in 
closed form and integrated to obtain 



For the Kolmogorov exponent m 2 = 0.773J)(rf). The 
intensity spectrum can also be used to show that the 
spatial scale of intensity in weak scattering is r{. As the 
strength of scintillation increases, refraction becomes im- 
portant and the intensity fluctuations show structure at 
two scales: a diffractive scale Sd.it, and a refractive scale 
s rol . In the limit of very strong scattering the electric 
field fluctuations become a zero-mean complex Gaussian 
random process, so the covariance of intensity becomes 
Tx(s) = |Fe(s)| 2 . Thus the e _1 scale of intensity s<jif is 
the e -1 ' 2 scale of the electric field so- This is the diffrac- 
tive limit. In the regime of moderately strong scintilla- 
tion the intensity can be modeled as unit variance diffrac- 
tive fluctuations modulated by refractive fluctuations as 
shown below (Rickett et al. 1984; Coles et al. 1987). 



l(t) = (l + 51 D (t))(l + 51 R (t)) 
= l + 51 D + SI R + 51 D 51 R . 



(4) 



The product term has the spatial scale of the more 
rapidly varying component 51d- The refractive term 
has a larger e~ 4 scale approximately the size of the 
scattering disc, s re f = 9qz = rf/so- The diffractive 
term has unit variance, and we define the variance of 



the refractive term as 



"ref- 



The total variance 



1 + 2m 2 cf and the total variance at the diffractive scale 
An asymptotic expression for the total vari- 



"ref- 



is 1 

ance in strong scintillation has been derived (Prokhorov 
et al. 1975), m 2 = 1 + (2(" +1 )/7ra) sin(7ra/2)r(l + 
a/2) 2 r(4/a - l)D(r [ ) ( -- ( - 2 / a ^ 2 - a ^ . For the Kolmogorov 
exponent m 2 = 1 + 0.476L>(r f )~°- 4 . 

In strong scattering the diffractive fluctuations are 
quite frequency dependent. One could estimate the co- 
variance Ti(Af) =< &T(fo)£T(fo + Af) >, which is equal 
to |r^(Af)| 2 in the Gaussian limit, because an expression 
for rg(Af) is available. However large low frequency fluc- 
tuations in phase decorrelate T^Af) and these must be 
removed. The theory of this is also discussed by Codona 
et al. (1986) but most observers use only the asymptotic 
result that Tx(Af) is the Fourier transform of the auto- 
correlation of I{t). Using the quadratic structure func- 
tion model this would give r z (Af) = 1/(1 + (2?rAf t ) 2 ), 
so the half power width of Tx is 5v = l/(27rfo)- From 
these definitions one can derive a number of useful re- 
lations and scaling factors for Kolmogorov spectra, e.g. 
Srcf/sdif = (r-f/so) 2 = io/Sv, ml ex A 2 ' 83 , s ex A" 1 ' 2 , 
s re f oc A 2,2 , and to oc A 4,4 . 

4. SIMULATED MEASUREMENT OF DIFFRACTIVE 
SCALE S D i F AND BANDWIDTH Sv 

A typical pulsar observation is reduced to a dynamic 
spectrum of integrated pulse power vs frequency and 
time. The observed time interval is usually longer than 
the diffractive time, but shorter than the refractive time. 
When one computes the two dimensional autocorrelation 
of the dynamic spectrum from such an observation, one 
obtains a biased estimate because the refractive intensity 
variations are not adequately sampled. A practical ques- 
tion is, "how does the undersampled refractive scintilla- 
tion affect the estimate of the diffractive scintillation?" 
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Fig. 1. — Dynamic spectra for m| = 10: (a) the top panel is the full dispersive simulation; (b) the middle panel is an expanded section of 
the top panel; (c) the bottom panel is the same as (b) except the simulation is run without dispersion. In each case the abscissa is scaled 
with respect to the Fresnel scale rf, (about 3.2 X 10 km for a typical pulsar). If converted to a time axis using a velocity of 100 km/s, the 
full window corresponds to about a month of observation. 



Wc have addressed this problem using simulated dy- 
namic spectra which are much longer than a typical ob- 
servation. A 16384 x 2048 (x,y) window was simulated 
and the center slice y — yo of length 16384 was extracted. 
Then the simulation was repeated at 128 different fre- 
quencies, building up a dynamic spectrum of 16384 x 128 
(x,f) from the central slices of each window. An exam- 
ple is shown in Figure la. In this particular simulation 
m 2 , = 10 and the abscissa has been scaled by rf = 20 
samples. The spatial scales derived from the measured 
covariance functions of electric field and intensity are 
so = 0.20 rf, Sdif = 0.22 rf and s re f = 5.5 rf (estimation of 
spatial scales will be discussed in detail in the following 
section). The theoretical scales are 0.215 rf, 0.215 rf and 
4.65 rf respectively. Thus s^n/dx = 4.4 so the smaller 
scales are well-sampled, and L y /s rc { = 19 so the largest 
scales fit easily in the window. For a typical pulsar at 
a distance of 200 pc with a scattering screen at 100 pc 
observed at 1.4 GHz the Fresnel scale rf = 3.2 x 10 5 km. 
The drift time for rf at 100 km/s would be about 1 hr, 
so the simulated window of 16384 samples corresponds 
to about a month of continuous observations. The same 
sampling was used for < to 2 < 30 but for ml > 30 we 
used 32768 x 4096 x 256. 

The expanded window plotted in Figure 1(b) shows 
both the diffractive and refractive structures more 
clearly. The expanded view also shows that the struc- 
ture is often tilted or even curved. This phenomenon, 



which will be discussed later, is primarily due to disper- 
sion. This is demonstrated by the bottom panel, Figure 
1(c), which shows the dynamic spectrum of the same 
phase screen without dispersion. 

Observers estimate the temporal and frequency scales 
of the observations by calculating the autocovariance of 
the dynamic spectrum. Of course this averages the tem- 
poral scale over the frequency range in the dynamic spec- 
trum, but the temporal scale varies almost linearly with 
frequency and the fractional bandwidth is usually less 
than 25% so this is not a significant bias. We computed 
the two dimensional autocovariance of the entire dynamic 
spectrum shown in Figure 1 (a). Cuts through this au- 
tocovariance are shown as solid lines on Figure 2. In the 
left panel Tx(s) is shown with a logarithmic abscissa so 
that one may see both the spatial scales clearly. The 
frequency correlation Ti(Af) is shown in the right panel. 
The fractional bandwidth 5v/io = 0.027, considerably 
smaller than the (so/ff) 2 = 0.046 expected theoretically. 
The diffractive scale is well-estimated because there are 
(ix/^dif) * (B/Sis) = 1.7 x 10 4 degrees of freedom in 
the dynamic spectrum (here B is the width of the dy- 
namic spectrum in frequency). However the refractive 
scale is broadband and much larger, so there are only 
(L x /s re {) * 1 = 136 degrees of freedom in the dynamic 
spectrum, and the refractive component is much less sta- 
ble. To improve our estimate of the refractive process we 
have computed the autocovariance of the entire 16384 x 



2048 window at the center frequency fo. This provides 
(Lx/sdii) * (Ly/sdii) — 2330 degrees of freedom and a 
much better estimate of the refractive component. 




Fig. 2. — Estimates of Vx(s) (left panel) and Tx(Af) (right panel) 
from the dynamic spectrum displayed in Figure 1(a) are shown as 
solid lines. In the left panel an estimate of Tx(s) calculated from 
the entire simulation plane at fo is also shown as a dashed line. It 
provides a more stable estimate of Tx(s) than does the dynamic 
spectrum alone, because there are more refractive scales in the 
entire simulation plane. The widths s^if and 5v/Iq are marked 
as solid circles. The abscissa in this panel is displayed with log 
scaling so both diffractive and refractive scales can be observed 
clearly. The zero points, which are not on the log abscissa, are 
marked as diamonds on the left axis. 

We have broken the full dynamic spectrum into 128 
blocks, each of length approximately s rc f, which corre- 
sponds to about 6.4 hrs of observation, and estimated the 
autocovariance in time and in frequency for each block. 
In each block there are 215 degrees of freedom in the 
diffractive process. From these we derived the spatial 
(temporal) and frequency scales for each block, by fitting 
a theoretical model to the appropriate autocovariance of 
each block. The results are shown in Figure 3. The error 
bars shown on Figure 3 are ±2<r as determined from the 
least squares fitting process. 

The theoretical model for the spatial correlation is 
Fx(s) = exp(— (s/sdif) 5 ' 3 ), which is the known asymp- 
totic form in strong scattering. The asymptotic forms 
for Fi(Af) are similar in shape but do not have a simple 
closed form (Coles et al., 1995b). We found that an expo- 
nential Ti(Af) = exp(— | ln(2)Af/<5z/|) gave a reasonable 
match to the theory and also to the overall average, so 
we used this model in the block fits. A gaussian model 
would give a significantly different fit. 

The spatial scale is shown in the upper panel of Fig- 
ure 3. The overall Sdif normalized to sq is 1.02, the 
weighted mean of the 128 blocks, normalized the same 
way, is 0.89. The unweighted rms of the 128 blocks is 
47% but the standard deviation computed from the er- 
ror bars estimated by the fitting process is only 9%, i.e. 
these error bars underestimate the actual error by a fac- 
tor of five. The frequency scale 8v normalized to fo is 
shown in the lower panel of Figure 3. The overall 8v /{$ 
is 0.027, whereas the weighted mean of the 128 blocks is 
0.017 so the mean of the blocks is only 63% of the over- 
all average. The unweighted rms of the blocks is 84% of 
the weighted mean but the standard deviation computed 
from the error bars is only 6%. Again the error bars se- 
riously underestimate the actual variation of the block 
estimates. 

Since the analysis shown in Figure 3 is typical of that 
used by observers, it is disappointing to find that the er- 
ror bars determined by a least squares fit are so unreliable 
in both time (space) scale and bandwidth. Observers 




Fig. 3. — Upper panel, apparent diffractive scale Sdif measured on 
blocks of length s rc f. Lower panel, apparent diffractive bandwidth 
Sv measured on the same blocks. The points marked with open 
circles are discussed below as examples of "good" and "bad" fits. 



have noted this effect and Cordes (1986) has suggested 
that the error bars should be computed with the num- 
ber of degrees of freedom reduced by a factor of 25 to 
100. In the simulation shown reducing the number of 
degrees of freedom from 215 to 2 would indeed match 
the estimated error to the observed rms variation. How- 
ever the variations do not have a Gaussian distribution, 
as demonstrated by the obvious "spikes". So there is a 
finite possibility of a much larger error. 

The spikes in time scale and bandwidth are not cor- 
related. The cause of the spikes is not obvious on an 
inspection of the block correlations. In Figure 4 we show 
six examples of the time scale fits. The top three panels 
are from the spikes at blocks 7, 85, and 126 which are 
marked with large open circles in Figure 3. The lower 
three panels are normal fits at samples 46, 47 and 48 
which are marked with small open circles in Figure 3. 
Although the top three panels represent huge errors, this 
is not apparent in the fit. Since this is a very impor- 
tant issue for observers, we have reanalyzed the data in 
blocks of half and twice the refractive scale. The results 
are provided in Table 1. As expected, the reported errors 
decrease with longer block lengths. However the rms is 
rather stable. This suggests that the rms is not caused by 
statistical errors in the fit, but by the refractive fluctua- 
tions. The downwards bias of the mean becomes worse 
with shorter blocks. Perhaps more important, the fre- 
quency of spikes is roughly independent of block length. 
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Fig. 4. — Least squares fit to correlation functions from six 
blocks. The top three panels are from spikes at blocks 7, 85, and 
126. The lower three panels are from more typical blocks 46, 47 and 
48. The noisy lines are the estimated covariances and the smooth 
lines are the best fit theoretical models. 



TABLE 1 

DlFFRACTIVE SCALE (NORMALIZED BY So) 



fit: block 


wtd mean 


rms 


error bar 


cov 2s re f 


0.99 


43% 


7% 


^rcf 


0.89 


47% 


9% 


Sref/2 


0.78 


53% 


12% 


spec s rc f 


0.72 


18% 


3% 



It is well known that least squares fits to correlation 
functions violate the normal assumptions that the er- 
rors (i.e. the deviations of the observations from the 
model) are independent and have equal variance. So we 
re-estimated the spatial scale by doing a least squares fit 
to the Fourier transform of the correlation function. The 
advantage of this is that the errors in the Fourier trans- 
form are independent and their variance is known, so it 
can be corrected in a weighted fit. We used exactly the 
same correlation model, but fit its Fourier transform to 
the Fourier transform of the correlation. The resulting 
spatial (time) scale is shown in Figure 5. 




Fig. 5. — Diffractive spatial (temporal) scale obtained by a 
weighted least squares fit in the frequency domain. 



One can see by eye that the spikes have completely dis- 
appeared, the rms is considerably smaller, but the mean 
is also lower. These parameters are also given in Table 
1. The weighted mean has decreased to 72% of the true 
value. The rms has decreased by a factor of three, but 
the error bars are still far below the rms. 

The same fits shown in Figure 4 are displayed in the 
frequency domain in Figure 6. One can see that there 
is a low frequency excess in the top panels, which gave 
rise to spikes when fit in the time domain, but these do 
not significantly distort the fit in the frequency domain 
because the optimal weighting favors the higher frequen- 
cies. Thus fitting the power spectrum is strongly rec- 
ommended. It provides a much better estimator of the 
spatial (temporal) scale, although users will have to cor- 
rect for a downwards bias of the order of 30% and the 
error bars will not be more reliable. Fitting in the spec- 
tral domain would be particularly valuable when one is 
analyzing the effect of orbital motion on scintillation of 
a binary pulsar, or when one is attempting to measure 
the annual modulation in time scale due to the Earth's 
velocity. It is clear that in all cases observers should be 
extremely cautious about the errors estimated from data 
blocks whose length is less than the refractive scale. 

We have not attempted to fit the frequency scale 8u in 
the Fourier domain, but we expect that it would also be 
better fit in this way. The theoretical model for Tj(Af) 
is demonstrably inadequate because the value of 8v mea- 
sured from the simulation shown is only 59% of that 
expected theoretically. The theory applicable to mod- 
erately strong scattering has been worked out (Codona 




Frequency 

Fig. 6. — Weighted least squares fit to the block power spectra. 
The noisy lines are the power spectra and the smooth lines are the 
best fit theoretical model. The blocks are the same as those shown 
in Figure 4. 



et al., 1986), but it has not been applied to interstel- 
lar scintillation. This should be done, but is beyond the 
scope of this paper. Observers have also noted this effect 
and Gupta et al. (1994) have given a heuristic model 
which explains the reduction in bv in terms of refrac- 
tive modulation of the diffractive scintillation, combined 
with the plasma dispersion. This model agrees approx- 
imately with the simulations. The bias changes slowly 
with increasing strength of scattering as will be shown 
later. 

We do not understand why the fluctuations in spa- 
tial scale are not correlated with those in bandwidth. 
This effect has been considered by Gupta et al, (1994) 
who suggested that refractive variations that tilt the spa- 
tial structure should change the bandwidth but not the 
scale. This is clearly not true of our simulations, both 
the bandwidth and the spatial scale are highly variable. 
These variations have also been observed in an exten- 
sive monitoring of ISS from a set of 18 pulsars by Bhat 
et al. (1999). They estimated the frequency and time 
scales fitting the covariance functions as we have done, 
but with a Gaussian model. We believe that the variabil- 
ity is due to the fact that some 20% of the variance is 
caused by refractive effects which are ignored, but which 
significantly modulate the diffractive effects. The vari- 
ous heuristic models employed by observers are useful, 
but this problem would benefit from further theoretical 
work based on Codona et al. (1986) and further simula- 
tions. In particular it would be very useful to extrapolate 
our results to stronger scintillation where the bias should 
decrease. 

5. SECONDARY SPECTRA 

It has been noticed for decades (Ewing et al. 1970; 
Roberts & Abies 1982), that dynamic spectra often 
show "criss-cross" structures. These structures, which 
are easily visible in Figure 1(b), are responsible for the 
"parabolic arcs" discovered by Stinebring et al. (2001) 
and discussed by Cordes et al. (2006), which are seen 
in the "secondary spectrum." The secondary spectrum, 
which is the two-dimensional power spectrum of the dy- 
namic spectrum and has been called the delay-Doppler 
spectrum, is shown in Figure 7. Here the delay axis has 
been scaled to ns assuming the center frequency is 1.4 
GHz. This scaling will be continued throughout. The 
Doppler axis is scaled by r{/V c g. With the values used 
earlier the Doppler range would be ±3 mHz. An arc is 



visible in Figure 7, but it is quite symmetric and not 
sharply defined. Most observed arcs are more clearly de- 
fined and they are often quite asymmetric. It is known 
that the arcs would be sharper in weaker scintillation 
or if the scattering were anisotropic (Cordes et al. 2006) 
and this is a strong argument that interstellar turbulence 
is often anisotropic. However it has not been understood 
exactly why observations often show asymmetric arcs. It 
has been assumed that this is due, in some way, to re- 
fraction, but it has been unclear whether this requires 
a discrete refracting structure, or whether the refraction 
that occurs naturally due to the low frequency part of 
the turbulent spectrum is sufficient. To test this question 
we have broken the full dynamic spectrum into smaller 
blocks, typical of the observation duration (as was done 
in calculating the time and frequency scales in Figure 
3). The results, which are shown in Figure 8, show ex- 
actly the characteristics of the observations. The arcs are 
sometimes quite asymmetric and the orientation of the 
asymmetry reverses on a time scale considerably longer 
than the refractive time scale. Evidently this resolves 
the question of whether such asymmetric secondary spec- 
tra are due to discrete deterministic structures or to re- 
fractive components of a turbulent spectrum - they can 
be caused by a continuous turbulent spectrum, but the 
scales involved are larger than the refractive scale. 




Differential Delay (ns at 1 .4 GHz) 



Fig. 7. — The delay-Doppler spectrum computed from the dy- 
namic spectrum shown in Figure 1 (a). The brightness scale is dB 
= lOlogio (power). 

When the same phase screen is simulated without dis- 
persion the dynamic spectrum, which is plotted in Fig- 
ure 1(c), shows criss-cross structures but very little tilt. 
The average secondary spectrum is very similar to that 
shown in Figure 7, but the secondary spectra of shorter 
blocks, which are shown in Figure 9, do not show the time 
varying asymmetry characteristic of the dispersive sim- 
ulation shown in Figure 8. Thus there can be no doubt 
that the phenomenon of tilted bands in dynamic spec- 
tra is not caused solely by refraction, but by dispersion 
and refraction. With this insight it is relatively simple 
to describe the phenomenon analytically. The physical 
difference can be explained intuitively. Large scale re- 
fraction will tilt an entire angular spectrum regardless of 
dispersion. However it will tilt to the minimum phase de- 
lay position (due to Fermat's Principle). The secondary 
spectrum is caused by interference of scattered waves 




Fig. 8. — Delay-Doppler distributions computed from 25 blocks 
of the dynamic spectrum shown in Figure 1(a). Time increases 
from the top left, first to the right, then down, ending at bottom 
right. Each panel has the same axes as Figure 7. 



with different group delay and different Doppler shift. 
In a non-dispersive problem the group delay is equal to 
the phase delay so a refractive shift of the entire angu- 
lar spectrum has no effect on the secondary spectrum. 
When the medium is dispersive the phase delay is the 
negative of the group delay, so the position of minimum 
phase delay is not the position of minimum group delay. 
Indeed it is a local maximum of group delay. In this case 
a refractive tilt has a dramatic effect on the parabolic 
arc. A rough analysis is given below. 




Fig. 9. — Delay-Doppler distributions computed from 25 blocks 
of the non-dispersive dynamic spectrum shown in Figure 1(c). 



We assume (without loss of generality) a thin screen 
half way between a pulsar and the observer. The arcs are 
caused by the interference of a wave scattered at an angle 
89 in the direction of the velocity and an unscattered 
wave. The two waves interfere with differential Doppler 
shifts 8fd = 2VS9/X and differential delays 8rd — zS9 2 /c 
where z is the distance from the screen to the observer. 
It is useful to normalize the parameters with respect to 
the rms scattering angle 6>o, i.e. fdo = 2V9 /X and Tdo — 
z9q/c. Then Srd/rdo — (Sfd/fdo) 2 defines the parabolic 
arc in normalized form. 

Now we consider the effect of a phase gradient V(/>* 
in the screen. It will tilt the angular spectrum seen by 
the observer by an angle 9* = \7<j)*/2k. If the gradient 



is perpendicular to the velocity then it will have no ef- 
fect on the arc. Thus we consider only the gradient in 
the direction of the velocity. The group delay is given by 
T g {9) = z8 2 /c=Fz#V(/>*/(x> in nondispersive (top sign) and 
dispersive (bottom sign) cases. The differential Dopplcr 
between an unscattered wave arriving at angle 9* and a 
scattered wave arriving at angle 8* + 59 is unchanged by 
the phase gradient. The differential delay is more com- 
plex. In the non-dispersive case we have 5t<i = zS9 2 /c, 
i.e. the arcs are unaffected by a phase gradient, at least 
to first order. However in the dispersive case we have 
6Td = z(S9 2 + 49*89) /c. Putting this in normalized form 
we have 8r d JT m = (Sf d /f d0 + 29*/9 ) 2 - (29*/9 ) 2 . The 
apex of the arc is shifted in Dopplcr by (29*/9o)fdo and in 
delay by (29*/9o) 2 Tdo. Thus, if the phase gradient shifts 
the entire angular spectrum by an amount comparable 
with the rms scattering angle 9q, then the shift in the 
apex of the parabola should be detectable in a disper- 
sive medium. This analysis shows why tilted arcs are ob- 
served in dispersive cases and not in non-dispersive cases. 
It also agrees with earlier rough analyses of the slope of 
tilted structures in the dynamic spectrum (Shishov 1974; 
Hewish 1980; Gupta et al. 1994). However, it is difficult 
to use these expressions quantitatively to estimate the 
electron density gradients in the IISM because the arcs 
are not always very distinct (as in the case simulated). 
To do this one would need to model the entire secondary 
spectrum including the effects of gradients both parallel 
and perpendicular to the velocity. 

6. TIME OF ARRIVAL FLUCTUATIONS 

6.1. Pulse Shape 

The simulated electric field E(x, y, /) allows one to 
calculate the pulse I(x,y,t) — \Ff{E(x,y, f)}\ 2 at each 
pixel. Here F{ is the Fourier transform operating on co- 
ordinate f. Thus the simulation provides a direct cal- 
culation of the pulse shape, including the arrival time 
variations. In the simulation the mean electron density 
over the simulation window is zero, but at any given po- 
sition there is "dispersion delay" and it fluctuates with 
position. In addition the pulse shape changes on both the 
diffractive and refractive scales. This is shown in Figure 
10. Here the pulse power is shown on a log 10 scale over 
a 40 dB dynamic range. The phase screen delay at f is 
ovcrplotted on Figure 10 as a white line. The screen de- 
lay is taken as the group delay of the phase screen, which 
is the negative of the phase delay for a plasma. In the top 
panel one can see that the peak of the pulse power tracks 
the slow variation due to dispersion measure. In the bot- 
tom panel one can see the much finer scale variation due 
to diffractive and refractive scattering. 

The average pulse shape would be dominated by the 
dispersion measure variations. However the scattered 
pulses can be aligned with respect to the dispersion de- 
lay in the screen itself. With this alignment all the delay 
variations are due to propagation from the screen to the 
observer. The average pulse-shape, which is the expected 
quasi-exponential, is shown in Figure 11. Note that this 
figure has a log ordinate, so an exponential pulse would 
drop linearly with delay. In order to display the leading 
edge of the pulse with low sidelobes we used a Black- 
man window in the Fourier transform for this display. 
One can see that the pulse tail does not drop nearly as 
fast as the exponential, which would be characteristic of 



a Gaussian angular spectrum. This is because the ac- 
tual angular spectrum falls more slowly than a Gaussian 
at large angles. It is easily shown that at high angles 
B{9) ex 9- ( - a+2 \ so I(t) ex £-(«+ 2 )/ 2 . The asymptotic 
diffractive theory is overplotted as a dashed line. The 
theoretical curve was scaled up by a factor of 1.5 to al- 
low for rounding of the peak by the Blackman window. 
It is clear that the long term average, when corrected for 
dispersion measure fluctuations, agrees very well with the 
simple diffractive theory. 

The diffractive effect is not to spread each pulse into 
a quasi-exponential, rather to break the pulse into sub- 
pulses. This is shown in a very expanded view in Figure 
12. It is only the superposition of all the pulses that is 
a continuous quasi-exponential. The width of the frag- 
ments of each pulse appears to be the resolution of the 
Fourier transform i.e. the inverse of the bandwidth. One 
never observes such breakup of the pulse shape in pul- 
sars because the intrinsic pulse width is always much 
larger than the inverse of the bandwidth. It might be 
observable in giant pulses with coherent de-dispersion 
(e.g. Hankins et al., 2003). 

6.2. Pulse Centroid 

Diffractive scattering also causes the centroid of the 
pulse to shift, and this effect is observable. Indeed it may 
be an important source of timing noise in some pulsars 
(Foster & Cordes 1990). The centroid is shown in the top 
panel of Figure 13 with the phase screen delay marked 
as a black line. One can see that the centroid of the 
pulse does not follow the screen delay as well in regions 
of high gradient in screen delay. This difference is much 
weaker in the non-dispersive case, which is shown in the 
lower panel. Here, of course, the screen delay plotted 
is the phase delay and it has the opposite sign of the 
group delay in the upper panel. Clearly a steep gradient 
in either case leads to increased refraction but the effect 
of this refraction on the delay is much greater in the 
dispersive case. 

It has been known for some time that scattering causes 
fluctuation in pulse arrival times, and that this fluctua- 
tion is anticorrelated with pulse power. This has been 
discussed theoretically (Blandford & Narayan 1985) and 
observed (Lestrade et al. 1998). Theory and observa- 
tions were discussed in terms of a refractive mechanism 
applied to observations which were averaged over many 
diffractive time scales. The anti-correlation also exists at 
diffractive scales, as can be seen in Figure 14 panels a and 
b. Here the first 50 rf of the simulation shown in Figure 
13 have been expanded to show the diffractive structure. 
The centroid corrected for the screen delay T c is shown 
in the top panel. In the middle panel we show the to- 
tal pulse flux density Pt ■ The anti-correlation is evident 
both at Sdif = 0.22rf and at s ro f = 5.5rf. The correlation 
is significantly higher between T c and l/y/Pr as appar- 
ent in Figure 14c. Here the raw cross correlation is 73%, 
and it rises to 83% if both series are low pass filtered to 
remove all scales larger than the refractive scale. 

The anti-correlation between TOA and flux at the 
diffractive scale has not been discussed theoretically, but 
one can see the mechanism involved by examining the 
shape of successive pulses on the diffractive time scale as 
shown in Figure 12. Evidently the pulse shape changes 
significantly on the diffractive time scale, so the instanta- 
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Fig. 10. — Pulse shape vs position. The intensity is loglO(powcr). The lower panel is an expanded section of the upper one. The white 
line is the group delay through the screen at the same transverse location as the intensity "measurement" . 
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Fig. 11. — Average pulse shape I(t), corrected for dispersion 
measure variations in the screen. The finite rise time is due to the 
bandwidth limit and the spectral window. The Blackman window 
is used because it provides very low sidelobes but it is considerably 
broader than a rectangular window. The dashed line is from the 
simple diffractivc theory. It has been shifted up by a factor of 1.5 
to allow for rounding of the peak by the Blackman window. 



neous angular spectrum must also vary significantly on 
that time scale. When the angular spectrum broadens 
the pulse broadens, and the flux drops correspondingly. 
However the pulse broadening is one-sided, it always in- 
creases the delay, so the delay of the centroid increases 
when the flux drops and we see this as an anti-correlation. 
One can see an interesting structure centered near 
x = 40rf, which very much resembles a structure ob- 




FlG. 12. — Expanded view of the pulse shape. Pulse delay runs 
downwards to the left. The total delay shown (80 ns) is comparable 
with the range shown on Figure 11. Observing time runs down- 
wards to the right. Pulse power is vertical. The period displayed 
is the drift time for 8 r{ or about 1.2 refractive time scales. 



served by Cognard et al. [f 993] who called it an extreme 
scattering event and attributed it to a discrete cloud of 
plasma. Our simulation shows that such events can occur 
naturally in Kolmogorov turbulence and this is confirmed 
by recent simulations very similar to ours [Hamidouche 
and Lestrade, 2007]. 

The T c fluctuations show a larger scale component (in 
addition to diffractive and refractive components) which 
is not correlated with \j\fPr- It does show a strong cor- 
relation with the magnitude of the gradient of the screen 
phase. This component is refractive but of larger scale 
than the refractive intensity fluctuations. Here the phase 
is essentially linear over the scattering disc, so we use a 
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Fig. 13. — Pulse ccntroid vs position. Top panel dispersive, 
bottom panel non-dispersive. The black line is the mean delay of 
the screen simulated at the observer's position. 




Distance x/r f 

Fig. 14. — The centroid corrected for dispersion measure T c (top 
panel). Total pulse power Pt (middle panel). 1/y/Pr (bottom 
panel). The cross correlation of the top and bottom figures is 73%. 
It rises to 83% when scales larger than 28 rf are filtered out. Here 
the first 50 rf of the 800 r f in the simulation have been expanded 
so the diffractive fluctuations can be seen clearly. 



plane wave approximation. A phase gradient V x 4> will 
cause an angular shift A9 X — V x <j)/k. This will displace 
the pattern by zA9 x = z\7 x (j)/k and thus cause a phase 
difference between the observed phase and the phase in 
the screen of A<fi — zA9 x V x (f> — (rfV x <t>) 2 ■ This causes 
a TOA error of A<fi/u) ~ (rfV x (f>) 2 /ui. There is also an 
incremental free space delay of half this TOA error, i.e. 
0.5 * A9 2 z/c — 0.5 * (rfV x (j)) 2 /uj. For a non-dispersive 
medium the free space delay partially cancels the refrac- 
tive component (an expression of Fermat's principle) and 
the cross correlation between T c and (rfV0) 2 /w would 
be negative. For dispersive medium the free space delay 
adds to the refractive component and the cross correla- 
tion is positive and larger. This crude analysis shows 
why a correlation between T c and (rfV0) 2 /w should be 
expected and, because the V y is ignored, why this cor- 
relation should not reach 100%. We will refer to this 
larger scale as the dispersive scale because it is stronger 
and positive in a dispersive simulation. The correlation 
can be seen clearly in Figure 15, where both T c and V0 
have been lowpass filtered to remove all the diffractive 
and refractive fluctuations. In this example the correla- 
tion coefficient for the dispersive case is 70% and for the 
non-dispersive case it is -37%. 

6.3. Correction of TOA Variation 

Since the TOA variations are highly correlated with in- 
tensity one can use the measured intensity to remove the 
correlated components of the TOA variation. We have 
done this in two steps. First we have reduced the correla- 
tion with l/y/Px to zero by subtracting a constant times 
1/y/Pr from T c , creating T cc . This reduced the variance 
in T c by more than a factor of two. We then reduced 
the correlation with (ffV 'if)) 2 / 'ui to zero in the same way, 
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Fig. 15. — A comparison of the corrected centroid from Figure 
14, smoothed over 28 rf to eliminate both refractive and diffractive 
fluctuations, with (rfV a; </') 2 /w where \7 x <t> has the same smooth- 
ing, (a) top panel the scattering medium is non dispersive; (b) 
middle panel the medium has the plasma dispersion; (c) the term 
('"fVicA) 2 /aj. The cross correlation of the middle and bottom pan- 
els is 70%. The cross correlation of the top and bottom panels is 
-37%. 



creating a T ccc . This reduced the variance by another 
20%. The three stages T c , T cc , and T ccc are shown in 
Figure 16 from top to bottom. Of course in practice our 
estimates of l/^/Pr and \d(DM)/dx\ will have error so 
the correction process will add some white noise, but in 
practice these errors are quite small compared with the 
white noise in T c so the added white noise is negligible. 
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Fig. 16. — The pulse centroid corrected for: (top panel) disper- 
sion measure fluctuations; (middle panel) intensity scintillation; 
(bottom panel) gradient of dispersion measure. The scales are the 
same in the three panels. One can see the effect of removing cor- 
related components clearly. 



7. CORRECTION OF TOA FLUCTUATIONS IN PRACTICE 

We have tested the correction process using observa- 
tions of J0437— 4715, a powerful nearby pulsar, at 685 
MHz, using the Parkes radio telescope with the CPSR2 
coherent de-disperser backend. The observations are sim- 
ilar to those described by Verbiest et al. (2008). These 
observations were made during a period of particularly 
heavy sampling between April 22, 2006 and June 11, 
2007. There were observations on 62 days with an aver- 
age of 140 minutes per day. The total IF bandwidth is 
64 MHz. Pulse arrival times were estimated on, roughly, 
5 minute subintegrations, and the signal to noise ratio 
(SN) was also estimated. The Verbiest et al. (2008) tim- 
ing model was used. The only timing parameter fitted 
was the pulse phase. Observations made simultaneously 
at 3 GHz were used to estimate the dispersion. The in- 
tegrated pulse power was not flux calibrated but it is 
believed that the 50 cm receiver noise is quite stable so 
the signal to noise ratio is proportional to the pulse flux. 
As this source is nearby the dispersion is relatively small, 
even at 685 MHz, but it is not negligible. The autocovari- 
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ance of 1/y/SN is shown in the left panel of Figure 17. 
It shows a quasi-exponential covariance with a time scale 
of about 20 min (which is the diffractive scale) and very 
little white noise (which would appear as a delta function 
at the origin). The autocovariance of the timing resid- 
uals corrected for dispersion is shown in the right panel 
of Figure 17. This covariance shows a clear white noise 
component and an exponential component with a time 
scale of about 30 min. This is not significantly different 
from the diffractive scale and is undoubtedly associated 
with the diffractive intensity variations. The white noise 
carries about the same variance as the 30 min variations. 
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Fig. 17. — Autocovariances of 1/%/SW (left panel) and TOA 
residuals (right panel) for pulsar J0437-4715 at 685 MHz. 



It should be noted that the rms white noise in these 
observations, about 1 /is, is much larger than the white 
noise in the normal observations of the Parkes Pulsar 
Timing Array (PPTA) for several reasons. The basic 
PPTA timing is done at 1400 MHz where the effects of 
the ISM are much smaller, with a much broader band- 
width and a much longer integration time, so the receiver 
noise is much lower. 

The cross correlation of the residuals with 1/y/SN is 
shown in Figure 18 before and after the scattering correc- 
tion. One can see that in advance of the scattering cor- 
rection the correlation (shown dotted) was about 50%. If 
one corrects for the white noise component in the residu- 
als the correlation of the diffractive component is about 
70% as in the simulations. After the scattering correc- 
tion the correlation (shown solid) is zero as expected. 
The scattering correction reduces the total variance of 
the residuals by about 25%. However one can see that 
the cross correlation persists as a negative component 
away from the origin. This could perhaps be removed 
using a more sophisticated removal. Rather than sub- 
tracting a constant times l/^/SN one could subtract the 
reference after being filtered with an FIR filter, as would 
be done in an adaptive echo canceller. 

It has been noted [Verbiest et al, 2009] that the errors 
in the residuals of J0437— 4715 do not behave like white 
noise. In particular when residuals sampled at 5 minute 
intervals are averaged the error on the average does not 
shrink like the square root of the number of residuals av- 
eraged. This is clearly because half the variance is carried 
in the diffractive component which has a correlation time 
considerably longer than 5 minutes. 

In this pulsar we are able to observe and correct diffrac- 
tive fluctuations in the TOAs. The refractive fluctuations 
were negligible but it seems likely that they too could 
have been corrected the same way. We attempted to an- 
alyze the refractive scattering contribution of the pulsar 
J1939+ 2134, as was done by Lestrade et al. [1998] how- 
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Fig. 18. — Cross correlation of TOA residuals with l/\/SN for 
pulsar J0437-4715 at 685 MHz. The trace shown dotted is before 
correction of the residuals for scattering and the trace shown solid 
is after this correction. 



ever the pulsar is not sampled adequately in the normal 
PPTA observations to resolve the 2.5 day time scale of 
the refractive scintillation. 

8. A HEURISTIC MODEL OF TOA VARIATION 

It would be very useful to be able to predict the TOA 
variations discussed above for a given pulsar at a given 
frequency. Since we lack an analytical theory, we have 
created a heuristic model with which we can scale the 
results of simulations to match the observing parameters. 
First we ran simulations over the range of strengths of 
scattering (to 2 ) accessible to our computing platforms: 
3, 5, 10, 20, 30, 60, and 100. Then we compared the 
behavior of the parameters for which we have theoretical 
scaling models, with the simulations. This comparison 
validates both the simulations and the scaling models. 

The fundamental spatial scales are: the e -05 scale of 
the electric field (sq); the diffractive e _1 scale of inten- 
sity (sdif); and the refractive e _1 scale of intensity (s re f). 
The scales of intensity were measured at their half power 
points, i.e. T = (3m 2 — l)/4 and (to 2 — l)/4, and then 
corrected to the e _1 values using the theoretical curve 
shape. The scale of the field is measured from the auto- 
covariance of the field at e~ ' 5 . The results are plotted 
vs to 2 in Figure 19 as symbols with error bars. The error 
bars are derived from multiple simulations. For a given 
m 2 we can use equation (3) to find sq. This expres- 
sion, which is valid in any strength of scattering, is plot- 
ted as a solid line and agrees well with the simulations. 
The strong scattering approximation for diffractive scale 
Sdif = so agrees less well but improves in stronger scatter- 
ing. The strong scattering approximation for refractive 
scale s rc f = r 2 /so, which is also plotted as a solid line, 
has weak agreement with the simulations. A heuristic 
expression for s 1G f derived from numerical solutions to 
the moment equations (Goodman and Narayan, 2006), 
which is plotted as dashed line, agrees much better. It 
is reassuring that the numerical solution agrees with the 
simulations for s rc f. 

The intensity variance is a fundamental parameter, but 
theoretical expressions are difficult to derive. We have 
plotted to 2 — 1 as symbols with error bars in Figure 
20. An asymptotic expression for large m 2 is plotted 
as a dashed line (Prokhorov et al. 1975). An empirical 
expression derived from numerical calculations is plot- 
ted as a solid line (Goodman and Narayan, 2006). One 
can see that the asymptotic expression converges very 
slowly, whereas the empirical expression is quite good 
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Born Variance mj; 
b 

Fig. 19. — Spatial scales from the simulations. The scales s^it 
(star), s re f (cross) and sn (circle) are normalized to r{. The error 
bars are derived from multiple simulations. The theoretical expres- 
sions drawn as solid lines fit so quite well, s^if less well, and s Ie { 
weakly. An empirical expression for s re f (Goodman and Narayan, 
2006) shown as a dashed line, fits well. 



throughout the simulated range. Observers should be 
cautious using the asymptotic expression because it evi- 
dently does not become accurate until m| > 100 and for 
such strong scattering the inner scale is likely to become 
important (increasing the variance) . Thus one should use 
the asymptotic approximation only when one is prepared 
to make a correction for the inner scale. 




Born Variance 



Fig. 20. — Intensity variance m? — 1 from the simulations. The 
error bars are derived from multiple simulations. The empirical 
expression (Goodman and Narayan, 2006) drawn as a solid line 
fits quite well. The asymptotic expression (Prokhorov et al., 1975) 
drawn as a dashed line converges very slowly. 



The fundamental time scale is the pulse width. A 
theoretical expression can be derived using so to obtain 
#o = 1/fcso and then in, — O^z/c. The measured pulse 
width (toobs) is obtained from the measured bandwidth 
8v by toobs = l/2irSv. The theoretical pulse width to and 
the measured TOA variations, normalized to the mea- 
sured pulse width ioo&s, are displayed in Figure 21. The 
ratio to/toobs is shown as square boxes. This ratio should 
be unity of course, but the dashed line through the data, 
which is given by to/toobs — 0.43(to^) 011 , is clearly a 
good approximation. We note that at ml = 10 the ob- 
served pulse width is twice the theoretical width. This 
was discussed in section 4 where we noted that the ob- 
served 8v is only half that theoretically predicted. This 
discrepancy decreases very slowly with increasing m 2 and 
probably becomes asymptotic to unity for m 2 >> 100. If 
this expression is used then toobs scales with A according 
to toobs oc A 409 . The wavelength dependence of toobs has 



long been a puzzle, as it has been found to vary more like 
A 4 , as is expected for a quadratic structure function, than 
the A 4 ' 4 which is predicted for a Kolmogorov spectrum. 
Our simulations show that this is simply a result of using 
the strong scattering formula outside its range of valid- 
ity. The simulations establish that ioobs can scale like 
A 41 for an isotropic Kolmogorov turbulence spectrum. 
Observations of this scaling behavior do not necessarily 
imply a steeper than Kolmogorov spectrum. 




Born Variance m b 



Fig. 21. — The theoretical pulse width and the rms centroid vari- 
ation for each simulation normalized to the observed pulse width 
toobs- The ratio ioAoobs i s marked with square boxes. The diffrac- 
tive, refractive, and dispersive contributions to the centroid vari- 
ation are marked with 'x', '+' and open circles respectively. The 
solid lines are empirical models with no theoretical justification. 



Estimates of the pulse arrival time are made by inte- 
grating over the duration of an observation T b s (typi- 
cally an hour for a PPTA observation) , and also over the 
observational bandwidth B (typically 256 MHz for a 20 
cm PPTA observation). The diffractive intensity fluctu- 
ations typically have a t&h < 1 h, and a 8v < 256 MHz. 
So one would expect some smoothing of the diffractive 
TOA fluctuations. Accordingly we have estimated the 
diffractive TOA variance with various bandwidths and 
we find that the rms diffractive TOA variation scales as 
(8v/ B) 1 / 2 as expected. We also find that the rms diffrac- 
tive TOA variation scales as (rdif/Tots) 1 ^ 2 . So in report- 
ing the results of the simulations we scale the diffractive 
TOA variations to B = Sv and T b s = rjif. These values 
can then be scaled using the actual observing time and 
bandwidth. The resulting values of T^/roobs are shown 
on Figure 21 as 'x' symbols. They are fit quite well with 
an empirical model Tdif = 3.0£oo6s(t*|) -0 ' 087 . 

The refractive time scale is generally much greater than 
the observing time, and the refractive fluctuations are 
correlated over the entire frequency band. Thus there is 
no significant smoothing of the refraction-induced TOA 
variations. Similarly the dispersive component of the 
TOA variations is not smoothed by the observational pa- 
rameters. These components T re f and T^is are shown on 
Figure 21 with '+' symbols and circles respectively. We 
do not have any theoretical basis for these parameters, 
but we expect them to be of the order of the pulse width 
because the pulse broadening is due to the superposition 
of these centroid variations. We extract empirical mod- 
els from these data for T re f = 1.62to bs(ml)~ ' 62 ', and 
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TABLE 2 
Predicted centroid variation for PPTA pulsars. The 

COLUMNS TO THE LEFT OF THE VERTICAL LINE ARE THE 
OBSERVATIONS ON WHICH THE PREDICTED CENTROID VARIATIONS 
TO THE RIGHT OF THE LINE ARE BASED. THE VALUES OF Su AND 

Tdif USED ARE THE GEOMETRIC MEAN OF THE OBSERVATION RANGE. 

THE PREDICTED VARIATION HAS BEEN CALCULATED AT 1400 MHz 
FOR THE FIRST 20 ROWS, AND AT 685 MHz FOR THE LAST ROW. 



source 


Jobs 


5v 


Tdif 


m b 


*0o6s 


Tdc 


T 
L ref 


x dis 


T 

^ rms 




MHz 


MHz 


min 




ns 


ns 


ns 


ns 


ns 


J0437-4715 


660 


17.0 


7.80 


1.94 


0.46 


0.74 


0.46 


0.23 


0.90 


J0613-0200 


1369 


1.74 


23.2 


188 


83.5 


5.60 


4.88 


2.01 


7.69 


J071 1-6830 


685 


2.32 


25.9 


11.7 


3.93 


2.95 


1.28 


0.59 


3.27 


J071 1-6830 


1369 


48.1 


67.6 


11.8 


3.03 


2.75 


0.98 


0.45 


2.96 


J1022+1001 


685 


16.0 


98.7 


2.34 


0.57 


2.79 


0.51 


0.25 


2.84 


J1024-0719 


685 


15.0 


45.2 


2.46 


0.61 


1.93 


0.52 


0.25 


2.02 


J1045-4509 


3100 


3.10 


5.14 


2324 


1234 


4.25 


15.1 


5.63 


16.7 


J1600-3053 


3100 


3.62 


9.59 


2041 


1056 


5.47 


14.0 


5.25 


16.0 


J1603-7202 


685 


1.68 


15.2 


15.3 


5.44 


2.55 


1.50 


0.68 


3.04 


J1603-7202 


1369 


5.69 


17.5 


70.0 


25.6 


3.12 


2.75 


1.18 


4.32 


J1643-1224 


3100 


1.33 


5.59 


4691 


2867 


6.07 


22.8 


8.23 


25.0 


J1713+0747 


685 


4.45 


30.7 


6.78 


2.05 


2.51 


0.94 


0.44 


2.72 


J1730-2304 


685 


2.55 


16.1 


10.8 


3.58 


2.25 


1.23 


0.57 


2.62 


J1730-2304 


1369 


11.2 


30.3 


39.9 


13.0 


3.18 


1.99 


0.87 


3.85 


J1732-5049 


1369 


3.34 


26.5 


109. 


43.6 


4.68 


3.56 


1.50 


6.07 


J1744-1144 


685 


12.8 


46.9 


2.81 


0.71 


2.09 


0.56 


0.27 


2.18 


J1824-2452 


3100 


0.81 


3.23 


7093 


4709 


5.57 


28.9 


10.3 


31.2 


J1857+0943 


685 


4.42 


16.9 


6.83 


2.07 


1.87 


0.94 


0.44 


2.14 


J1857+0943 


1369 


8.22 


33.0 


51.5 


17.7 


3.73 


2.31 


1.00 


4.50 


J1909-3744 


685 


7.48 


33.2 


4.40 


1.22 


2.15 


0.73 


0.35 


2.30 


J1939+2134 


1369 


2.20 


6.40 


154. 


66.0 


2.69 


4.36 


1.80 


5.43 


J2124-3358 


436 


6.90 


44.0 


0.90 


0.22 


1.74 


0.35 


0.18 


1.78 


J2129-5718 


685 


2.77 


21.6 


10.1 


3.30 


2.52 


1.18 


0.55 


2.84 


J2129-5718 


1369 


76.5 


52.6 


8.03 


1.90 


2.04 


0.78 


0.37 


2.21 


J2145-0750 


685 


11.3 


47.1 


3.12 


0.81 


2.19 


0.60 


0.29 


2.29 


J0437-4715 


660 


17.0 


7.80 


14.7 


8.07 


2.98 


2.29 


1.04 


2.62 



Tdis = 0.87£oo6s("i6)~ ' 66 - These curves are plotted over 
the simulated points as solid lines. The dispersive simu- 
lation for ml = 100 seems to be out of character with the 
rest of the data. This is probably because the simulated 
screen was not long enough to separate the refractive 
and dispersive components cleanly. As this screen was 
already 1 GB in size, we couldn't increase it easily. 

9. COMPARISON WITH PPTA OBSERVATIONS 

We have used the observed diffractive scintillation pa- 
rameters of the pulsars in the PPTA to estimate the TOA 
fluctuations due to scattering. For all the pulsars mea- 
surements of Tdif and 5v are available at some frequency 
foi> Sl although they are sometimes quite noisy and ap- 
parently variable. These parameters are seldom mea- 
sured at 1400 MHz, which is the primary frequency at 
which TOAs are measured with the PPTA, so we scale 
them all to 1400 MHz first. We use the Kolmogorov 
scaling for Tdif oc f ' , but we scale 5v oc r"° as dis- 
cussed in the previous section. From these we can cal- 
culate ml = 0.773(f/5^) 5 / 6 and t 0ohs = \/2tt5v. We 
then find T dif , T rcf , and T dis at 1400 MHz using the 
empirical models discussed in the previous section. Fi- 
nally we use Tdif and 8v to estimate the number of in- 
dependent diffractive scintles in the observation NS = 
(?obs/Tdif)(B/<5^). The diffractive contribution to the ob- 
served TOA is T dc = T di{ /VNS- The total predicted rms 
of the TOA is obtained by adding the three components 
in quadrature. The results are shown in Table 2. 



The predicted scattering noise at 1400 MHz is less than 
50 ns rms for all of the pulsars in the regular PPTA ob- 
servations. If this prediction is accurate then scattering 
is not a significant source of timing noise at the PPTA. 
However we know of two observations of timing noise 
which is correlated with intensity variations. The first we 
have already discussed, the timing noise in J0437— 4715 
at 685 MHz is correlated at the diffractive time scale. 
The observed rms is 1000 ns and the predicted rms is 
only 2.6 ns. In this case the expected pulse width is only 
8.07 ns, so it is difficult to understand how scintillation 
could cause an rms 100 times larger than the pulse width. 
The second case of such correlation has been observed in 
the pulsar J1939+2134 at the refractive time scale at 1.4 
GHz (Lestrade et al 1998). Here the rms TOA residuals 
were of the order of 300 ns and the autocorrelation was of 
the order of 25% at the first time lag. This is consistent 
with about 150 ns of timing noise correlated with flux. 
However we only predict 5.4 ns of timing noise, which 
would not have been observable. The predicted pulse 
width is only 66 ns, so it is surprising to see such a large 
observed variation. In these observations the diffractive 
scale could not be observed, so the correlated rms could 
be even higher. 

The simulations reported here are for a thin scattering 
layer of homogeneous isotropic turbulence with a Kol- 
mogorov spectrum. The two pulsars in question are well 
approximated by a thin layer of Kolmogorov turbulence 
(Ramachandran et al. 2006; You ct al. 2007), but the 
turbulence may not be isotropic or homogeneous. Indeed 
it has been suggested that J1939+2134 underwent an ex- 
treme scattering event (Cognard et al. 1993; Lestrade et 
al. 1998). It is becoming more evident that many pulsars 
show the effects of anisotropic scattering and it has been 
shown that anisotropy can have a pronounced effect on 
pulse broadening (Rickett et al. 2009), but we do not 
have an estimate of the anisotropy of the scattering for 
either J0437-4715 or J1939+2134. It has recently been 
shown that the scattering in a different pulsar B0834+06 
is caused by highly inhomogeneous turbulence (Brisken 
et al. 2009) and this may also be implied by the earlier 
observations of Hill et al. (2005). It is easy to imagine 
that inhomogeneous turbulence would greatly enhance 
refractive scattering but it is much less obvious what it 
might do to diffractive effects. The simulation can be 
extended to study the effects of anisotropy and inhomo- 
geneity and we are undertaking such an extension. 

It is also possible that there is an unrecognized mech- 
anism which causes timing fluctuations which are cor- 
related with signal to noise ratio. If so, it should be 
correctable in the same way, but it is very important to 
understand the mechanism. 

10. CONCLUSIONS 

Simulations have been shown to be a very useful tool 
in understanding observations of scattering in the inter- 
stellar plasma. The simulation engine has been regular- 
ized to the point where it can be used by others and it 
is available from the authors. Here we show a number 
of examples where simulations have allowed us to un- 
derstand long standing peculiarities in the observations. 
We also show that the simulated timing noise due to ho- 
mogeneous isotropic turbulence is considerably smaller 
than two cases of observed timing noise which is corre- 
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lated with intensity fluctuations. This opens a number 
of interesting questions about the effects of the assump- 
tions of homogeneous isotropic turbulence, and about the 
observations themselves. 
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